#############################################################
## Estimating Effects of English Rule on Litigation Outcomes
## Eric Helland and Jungmo Yoon (2016)
## This file has a few R functions
## Please read our web appendix for more information   

MH.bound <- function(y1,y0,ts,tau,never=0){
	# Estimating quantile treatment effects at tau
	# never=0 for trial outcomes, never=1 for settlemt outcomes
	# y1 (y0) is the outcome under the English (American) rule
	# ts is the trimming fraction
	if(never==0){
		QY1.U <- quantile(y1,(ts+(1-ts)*tau))
		QY1.L <- quantile(y1,((1-ts)*tau))
		QY0   <- quantile(y0,tau)
		QY1   <- c(QY1.L,QY1.U)
		Diff  <- QY1 - rep(QY0,2)
	}
	if(never==1){
		QY0.U <- quantile(y0,(ts+(1-ts)*tau))
		QY0.L <- quantile(y0,((1-ts)*tau))
		QY1   <- quantile(y1,tau)
		QY0   <- c(QY0.L,QY0.U)
		Diff  <- rep(QY1,2) - rev(QY0)
	}
	return(list(QY1 = QY1, QY0 = QY0, QTE = Diff))
}

Lee.bound <- function(y1,y0,ps,never=0){
	# Estimating the average treatment effect
	# never=0 for trial outcome, never=1 for settlemt outcome
	# y1 (y0) is the outcome under the English (American) rule
	# ps is the trimming fraction
	if(never==0){
		yp.L  <- quantile(y1,(1-ps))
		yp.U  <- quantile(y1,ps)
		AY1.L <- mean(y1[y1 <= yp.L])
		AY1.U <- mean(y1[y1 > yp.U])
		AY0   <- mean(y0)
		AY1   <- c(AY1.L,AY1.U)
		Diff  <- AY1 - rep(AY0,2)
	}
	if(never==1){
		yp.L  <- quantile(y0,(1-ps))
		yp.U  <- quantile(y0,ps)
		AY0.L <- mean(y0[y0 <= yp.L])
		AY0.U <- mean(y0[y0 > yp.U])
		AY1   <- mean(y1)
		AY0   <- c(AY0.L,AY0.U)
		Diff  <- rep(AY1,2) - rev(AY0)
	}
	return(list(AY1 = AY1, AY0 = AY0, ATE = Diff))
}

MH.bound.se <- function(y1,y0,ts,nn,a0,E.d,E.sd,E.s1d,E.1sd,E.1s1d,tau,never=0){
	# Standard errors for quantile effects 
	# never=0 for trial outcomes, never=1 for settlemt outcomes	
	if(never==0){
		QY1.U <- quantile(y1,(ts+(1-ts)*tau))
		QY1.L <- quantile(y1,((1-ts)*tau))
		QY0   <- quantile(y0,tau)
		QY1   <- c(QY1.L,QY1.U)
		Diff  <- QY1 - rep(QY0,2)	
		tau.l <- tau*(1-ts)
		tau.u <- tau*(1-ts)+ts
		W.c <- (tau*(1-tau))/((eval.f(y0,QY0)^2)*E.s1d)
		W.p <- ((a0/(1-ts))*(1-(a0/(1-ts))))/(((a0/((1-ts)^2))^2)*E.d) + (a0*(1-a0))/(((a0/(1-ts))^2)*(1-E.d))
		W.lb1 <- (tau.l*(1-tau.l))/((eval.f(y1,QY1.L)^2)*E.sd)
		W.lb2 <- (tau^2)/(eval.f(y1,QY1.L)^2)*W.p
		W.lb  <- (1/nn)*(W.lb1 + W.lb2 + W.c)
		W.ub1 <- (tau.u*(1-tau.u))/((eval.f(y1,QY1.U)^2)*E.sd)
		W.ub2 <- ((1-tau)^2)/(eval.f(y1,QY1.U)^2)*W.p
		W.ub  <- (1/nn)*(W.ub1 + W.ub2 + W.c)
		se.bnd<- c(sqrt(W.lb),sqrt(W.ub))
	}
	if(never==1){
		QY0.U <- quantile(y0,(ts+(1-ts)*tau))
		QY0.L <- quantile(y0,((1-ts)*tau))
		QY1   <- quantile(y1,tau)
		QY0   <- c(QY0.L,QY0.U)
		Diff  <- rep(QY1,2) - rev(QY0)	
		tau.u <- tau*(1-ts)
		tau.l <- tau*(1-ts)+ts
		W.c <- (tau*(1-tau))/((eval.f(y1,QY1)^2)*E.1sd)
		W.p <- ((a0/(1-ts))*(1-(a0/(1-ts))))/(((a0/((1-ts)^2))^2)*(1-E.d)) + (a0*(1-a0))/(((a0/(1-ts))^2)*E.d)
		W.lb1 <- (tau.l*(1-tau.l))/((eval.f(y0,QY0[1])^2)*E.1s1d)
		W.lb2 <- ((1-tau)^2)/(eval.f(y0,QY0[1])^2)*W.p
		W.lb  <- (1/nn)*(W.lb1 + W.lb2 + W.c)
		W.ub1 <- (tau.u*(1-tau.u))/((eval.f(y0,QY0[2])^2)*E.1s1d)
		W.ub2 <- (tau^2)/(eval.f(y0,QY0[2])^2)*W.p
		W.ub  <- (1/nn)*(W.ub1 + W.ub2 + W.c)
		se.bnd<- c(sqrt(W.lb),sqrt(W.ub))
	}	
	return(list(QY1 = QY1, QY0 = QY0, QTE = Diff, se.QTE = se.bnd))
}

Lee.bound.se <- function(y1,y0,ps,nn,a0,E.d,E.sd,E.s1d,E.1sd,E.1s1d,never=0){
	# Estimating the average treatment effect
	# never=0 for trial outcome, never=1 for settlemt outcome
	if(never==0){
		yp.L  <- quantile(y1,(1-ps))
		yp.U  <- quantile(y1,ps)
		AY1.L <- mean(y1[y1 <= yp.L])
		AY1.U <- mean(y1[y1 > yp.U])
		AY0   <- mean(y0)
		AY1   <- c(AY1.L,AY1.U)
		Diff  <- AY1 - rep(AY0,2)
		var.L <- var(y1[y1 <= yp.L])
		var.U <- var(y1[y1 > yp.U])
		V.c <- var(y0)/E.s1d
		V.p <- ((1-ps)^2)*((1-(a0/(1-ps)))/(E.d*(a0/(1-ps)))+(1-a0)/((1-E.d)*a0))
		V.lb1 <- var.L/(E.sd*(1-ps)) + ((yp.L-AY1.L)^2)*ps/(E.sd*(1-ps)) + (((yp.L-AY1.L)/(1-ps))^2)*V.p
		V.ub1 <- var.U/(E.sd*(1-ps)) + ((yp.U-AY1.U)^2)*ps/(E.sd*(1-ps)) + (((yp.U-AY1.U)/(1-ps))^2)*V.p
		V.lb  <- (1/nn)*(V.lb1 + V.c)
		V.ub  <- (1/nn)*(V.ub1 + V.c)
		se.bnd<- c(sqrt(V.lb),sqrt(V.ub))
	}
	if(never==1){
		yp.L  <- quantile(y0,ps)
		yp.U  <- quantile(y0,(1-ps))
		AY0.L <- mean(y0[y0 > yp.L])
		AY0.U <- mean(y0[y0 <= yp.U])
		AY1   <- mean(y1)
		AY0   <- c(AY0.L,AY0.U)
		Diff  <- rep(AY1,2) - AY0
		var.L <- var(y0[y0 > yp.L])
		var.U <- var(y0[y0 <= yp.U])
		V.c <- var(y1)/E.1sd
		V.p <- ((1-ps)^2)*((1-(a0/(1-ps)))/((1-E.d)*(a0/(1-ps)))+(1-a0)/(E.d*a0))
		V.lb1 <- var.L/(E.1s1d*(1-ps)) + ((yp.L-AY0.L)^2)*ps/(E.1s1d*(1-ps)) + (((yp.L-AY0.L)/(1-ps))^2)*V.p
		V.ub1 <- var.U/(E.1s1d*(1-ps)) + ((yp.U-AY0.U)^2)*ps/(E.1s1d*(1-ps)) + (((yp.U-AY0.U)/(1-ps))^2)*V.p
		V.lb  <- (1/nn)*(V.lb1 + V.c)
		V.ub  <- (1/nn)*(V.ub1 + V.c)
		se.bnd<- c(sqrt(V.lb),sqrt(V.ub))
	}
	return(list(AY1 = AY1, AY0 = rev(AY0), ATE = Diff, se.ATE = se.bnd))
}

eval.f <- function(ye,qe){
	# used to calculate the standard errors
	f1 <- density(ye,adjust=2)
	fq <- approx(f1$x,f1$y,xout=qe,method="linear")$y
	return(fq)
}
